Load libraries

Read in the data first, then combine Woodcock and Wiat and z-score the Standardized Test composite

d <- read.csv("/Users/jsulliv1/Downloads/jesstimation-master/zenith all data.csv")
d$condition <- factor(d$condition, levels=c("control","abacus"))
d$standardized <- rowMeans(d[,c("wiat","woodcock")],na.rm=TRUE) 

Descriptives and preliminaries

Summaries.

d %>% 
  select(standardized, ans, deviance, linr2, ordinality, year) %>%
  group_by(year) %>% 
  do(tidy(summary(.))) %>% 
  data.frame
##     year Var1          Var2              Freq
## 1      0       standardized Min.   :0.04651  
## 2      0       standardized 1st Qu.:0.15302  
## 3      0       standardized Median :0.19302  
## 4      0       standardized Mean   :0.19596  
## 5      0       standardized 3rd Qu.:0.23791  
## 6      0       standardized Max.   :0.34279  
## 7      0       standardized              <NA>
## 8      0                ans  Min.   :0.1000  
## 9      0                ans  1st Qu.:0.2325  
## 10     0                ans  Median :0.3600  
## 11     0                ans  Mean   :0.3657  
## 12     0                ans  3rd Qu.:0.4600  
## 13     0                ans  Max.   :0.7800  
## 14     0                ans      NA's   :25  
## 15     0           deviance     Min.   : NA  
## 16     0           deviance     1st Qu.: NA  
## 17     0           deviance     Median : NA  
## 18     0           deviance     Mean   :NaN  
## 19     0           deviance     3rd Qu.: NA  
## 20     0           deviance     Max.   : NA  
## 21     0           deviance     NA's   :187  
## 22     0              linr2     Min.   : NA  
## 23     0              linr2     1st Qu.: NA  
## 24     0              linr2     Median : NA  
## 25     0              linr2     Mean   :NaN  
## 26     0              linr2     3rd Qu.: NA  
## 27     0              linr2     Max.   : NA  
## 28     0              linr2     NA's   :187  
## 29     0         ordinality     Min.   : NA  
## 30     0         ordinality     1st Qu.: NA  
## 31     0         ordinality     Median : NA  
## 32     0         ordinality     Mean   :NaN  
## 33     0         ordinality     3rd Qu.: NA  
## 34     0         ordinality     Max.   : NA  
## 35     0         ordinality     NA's   :187  
## 36     0               year       Min.   :0  
## 37     0               year       1st Qu.:0  
## 38     0               year       Median :0  
## 39     0               year       Mean   :0  
## 40     0               year       3rd Qu.:0  
## 41     0               year       Max.   :0  
## 42     0               year              <NA>
## 43     1       standardized Min.   :0.07814  
## 44     1       standardized 1st Qu.:0.26279  
## 45     1       standardized Median :0.31116  
## 46     1       standardized Mean   :0.31115  
## 47     1       standardized 3rd Qu.:0.36023  
## 48     1       standardized Max.   :0.54093  
## 49     1       standardized              <NA>
## 50     1                ans  Min.   :0.0200  
## 51     1                ans  1st Qu.:0.1200  
## 52     1                ans  Median :0.1500  
## 53     1                ans  Mean   :0.1773  
## 54     1                ans  3rd Qu.:0.2000  
## 55     1                ans  Max.   :0.7800  
## 56     1                ans       NA's   :6  
## 57     1           deviance   Min.   :0.230  
## 58     1           deviance   1st Qu.:0.480  
## 59     1           deviance   Median :0.605  
## 60     1           deviance   Mean   :0.708  
## 61     1           deviance   3rd Qu.:0.880  
## 62     1           deviance   Max.   :2.400  
## 63     1           deviance       NA's   :3  
## 64     1              linr2  Min.   :0.0100  
## 65     1              linr2  1st Qu.:0.2300  
## 66     1              linr2  Median :0.3450  
## 67     1              linr2  Mean   :0.3663  
## 68     1              linr2  3rd Qu.:0.5100  
## 69     1              linr2  Max.   :0.9100  
## 70     1              linr2       NA's   :3  
## 71     1         ordinality  Min.   :0.5900  
## 72     1         ordinality  1st Qu.:0.7575  
## 73     1         ordinality  Median :0.8000  
## 74     1         ordinality  Mean   :0.7977  
## 75     1         ordinality  3rd Qu.:0.8400  
## 76     1         ordinality  Max.   :1.0000  
## 77     1         ordinality       NA's   :3  
## 78     1               year       Min.   :1  
## 79     1               year       1st Qu.:1  
## 80     1               year       Median :1  
## 81     1               year       Mean   :1  
## 82     1               year       3rd Qu.:1  
## 83     1               year       Max.   :1  
## 84     1               year              <NA>
## 85     2       standardized  Min.   :0.1298  
## 86     2       standardized  1st Qu.:0.3744  
## 87     2       standardized  Median :0.4251  
## 88     2       standardized  Mean   :0.4258  
## 89     2       standardized  3rd Qu.:0.4801  
## 90     2       standardized  Max.   :0.6358  
## 91     2       standardized       NA's   :1  
## 92     2                ans  Min.   :0.0100  
## 93     2                ans  1st Qu.:0.1100  
## 94     2                ans  Median :0.1300  
## 95     2                ans  Mean   :0.1462  
## 96     2                ans  3rd Qu.:0.1700  
## 97     2                ans  Max.   :0.5700  
## 98     2                ans       NA's   :2  
## 99     2           deviance  Min.   :0.2800  
## 100    2           deviance  1st Qu.:0.4500  
## 101    2           deviance  Median :0.6100  
## 102    2           deviance  Mean   :0.7127  
## 103    2           deviance  3rd Qu.:0.8850  
## 104    2           deviance  Max.   :3.0700  
## 105    2           deviance              <NA>
## 106    2              linr2  Min.   :0.0600  
## 107    2              linr2  1st Qu.:0.2050  
## 108    2              linr2  Median :0.3200  
## 109    2              linr2  Mean   :0.3518  
## 110    2              linr2  3rd Qu.:0.4700  
## 111    2              linr2  Max.   :0.8000  
## 112    2              linr2              <NA>
## 113    2         ordinality  Min.   :0.6200  
## 114    2         ordinality  1st Qu.:0.7500  
## 115    2         ordinality  Median :0.8000  
## 116    2         ordinality  Mean   :0.7936  
## 117    2         ordinality  3rd Qu.:0.8400  
## 118    2         ordinality  Max.   :0.9700  
## 119    2         ordinality              <NA>
## 120    2               year       Min.   :2  
## 121    2               year       1st Qu.:2  
## 122    2               year       Median :2  
## 123    2               year       Mean   :2  
## 124    2               year       3rd Qu.:2  
## 125    2               year       Max.   :2  
## 126    2               year              <NA>
## 127    3       standardized  Min.   :0.1614  
## 128    3       standardized  1st Qu.:0.4609  
## 129    3       standardized  Median :0.5442  
## 130    3       standardized  Mean   :0.5381  
## 131    3       standardized  3rd Qu.:0.6216  
## 132    3       standardized  Max.   :0.7740  
## 133    3       standardized              <NA>
## 134    3                ans  Min.   :0.0600  
## 135    3                ans  1st Qu.:0.1000  
## 136    3                ans  Median :0.1200  
## 137    3                ans  Mean   :0.1384  
## 138    3                ans  3rd Qu.:0.1600  
## 139    3                ans  Max.   :0.5900  
## 140    3                ans       NA's   :2  
## 141    3           deviance  Min.   :0.2300  
## 142    3           deviance  1st Qu.:0.4300  
## 143    3           deviance  Median :0.6250  
## 144    3           deviance  Mean   :0.6984  
## 145    3           deviance  3rd Qu.:0.9000  
## 146    3           deviance  Max.   :2.0200  
## 147    3           deviance       NA's   :1  
## 148    3              linr2  Min.   :0.0400  
## 149    3              linr2  1st Qu.:0.1925  
## 150    3              linr2  Median :0.3100  
## 151    3              linr2  Mean   :0.3560  
## 152    3              linr2  3rd Qu.:0.4975  
## 153    3              linr2  Max.   :0.8200  
## 154    3              linr2       NA's   :1  
## 155    3         ordinality  Min.   :0.5700  
## 156    3         ordinality  1st Qu.:0.7525  
## 157    3         ordinality  Median :0.8000  
## 158    3         ordinality  Mean   :0.7982  
## 159    3         ordinality  3rd Qu.:0.8500  
## 160    3         ordinality  Max.   :0.9500  
## 161    3         ordinality       NA's   :1  
## 162    3               year       Min.   :3  
## 163    3               year       1st Qu.:3  
## 164    3               year       Median :3  
## 165    3               year       Mean   :3  
## 166    3               year       3rd Qu.:3  
## 167    3               year       Max.   :3  
## 168    3               year              <NA>

Reliability for Estimation tasks

Deviance.

data.frame(year=factor(c("1-2","2-3")),
           corrs=c(cor.test(d$deviance[d$year==1],
                            d$deviance[d$year==2])$estimate,
                   cor.test(d$deviance[d$year==2],
                            d$deviance[d$year==3])$estimate))
##   year     corrs
## 1  1-2 0.2775883
## 2  2-3 0.5415512
data.frame(year=factor(c("1-2","2-3")),
           pvals=c(cor.test(d$deviance[d$year==1],
                            d$deviance[d$year==2])$p.value,
                   cor.test(d$deviance[d$year==2],
                            d$deviance[d$year==3])$p.value))
##   year        pvals
## 1  1-2 1.361945e-04
## 2  2-3 1.332268e-15

Linear \(r^2\).

data.frame(year=factor(c("1-2","2-3")),
           corrs=c(cor.test(d$linr2[d$year==1],d$linr2[d$year==2])$estimate,
                   cor.test(d$linr2[d$year==2],d$linr2[d$year==3])$estimate))
##   year     corrs
## 1  1-2 0.1907387
## 2  2-3 0.3790020
data.frame(year=factor(c("1-2","2-3")),
          pvals=c(cor.test(d$linr2[d$year==1],d$linr2[d$year==2])$p.value,
                   cor.test(d$linr2[d$year==2],d$linr2[d$year==3])$p.value))
##   year        pvals
## 1  1-2 9.499658e-03
## 2  2-3 9.585397e-08

Ordinality.

data.frame(year=factor(c("1-2","2-3")),
           corrs=c(cor.test(d$ordinality[d$year==1],
                            d$ordinality[d$year==2])$estimate,
                   cor.test(d$ordinality[d$year==2],
                            d$ordinality[d$year==3])$estimate))
##   year     corrs
## 1  1-2 0.2130936
## 2  2-3 0.3334274
data.frame(year=factor(c("1-2","2-3")),
          pvals=c(cor.test(d$ordinality[d$year==1],
                            d$ordinality[d$year==2])$p.value,
                   cor.test(d$ordinality[d$year==2],
                            d$ordinality[d$year==3])$p.value))
##   year        pvals
## 1  1-2 3.681035e-03
## 2  2-3 3.311927e-06

ANS.

data.frame(year=factor(c("0-1","1-2","2-3")),
           corrs=c(cor.test(d$ans[d$year==0],
                            d$ans[d$year==1])$estimate, 
                   cor.test(d$ans[d$year==1],
                            d$ans[d$year==2])$estimate,
                   cor.test(d$ans[d$year==2],
                            d$ans[d$year==3])$estimate))
##   year     corrs
## 1  0-1 0.2318435
## 2  1-2 0.1907536
## 3  2-3 0.4443862
data.frame(year=factor(c("0-1","1-2","2-3")),
           pvals=c(cor.test(d$ans[d$year==0],
                            d$ans[d$year==1])$p.value, 
                   cor.test(d$ans[d$year==1],
                            d$ans[d$year==2])$p.value,
                   cor.test(d$ans[d$year==2],
                            d$ans[d$year==3])$p.value))
##   year        pvals
## 1  0-1 3.377277e-03
## 2  1-2 1.010514e-02
## 3  2-3 2.630023e-10

We see that the correlations are pretty decent year-to-year for our three DVs.

Training effects

OK, now we test for effect of training for ANS. We can use a model b/c we have all 4 years:

w.interaction <- lmer(ans ~ year  * condition + (subnum|year), data=d)
summary(w.interaction)
## Linear mixed model fit by REML ['lmerMod']
## Formula: ans ~ year * condition + (subnum | year)
##    Data: d
## 
## REML criterion at convergence: -1205.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7297 -0.4898 -0.1617  0.3022  6.0360 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.  Corr 
##  year     (Intercept) 4.513e-03 6.718e-02      
##           subnum      1.367e-12 1.169e-06 -1.00
##  Residual             1.022e-02 1.011e-01      
## Number of obs: 713, groups:  year, 4
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)           0.324419   0.056844   5.707
## year                 -0.074453   0.030366  -2.452
## conditionabacus      -0.021833   0.013056  -1.672
## year:conditionabacus  0.006231   0.006859   0.908
## 
## Correlation of Fixed Effects:
##             (Intr) year   cndtnb
## year        -0.802              
## conditinbcs -0.113  0.091       
## yr:cndtnbcs  0.092 -0.110 -0.814
wo.interaction<-lmer(ans ~ year + condition + (subnum|year), data=d)
anova(w.interaction, wo.interaction)
## refitting model(s) with ML (instead of REML)
## Data: d
## Models:
## wo.interaction: ans ~ year + condition + (subnum | year)
## w.interaction: ans ~ year * condition + (subnum | year)
##                Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## wo.interaction  7 -1217.2 -1185.2 615.62  -1231.2                         
## w.interaction   8 -1216.0 -1179.5 616.03  -1232.0 0.8209      1     0.3649

Now, we do t-tests by task to find differences between control and abacus groups. We will correct for multiple comparisons post-hoc, correcting for the # of comparisons per DV.

tasks <- c("deviance","linr2","ans", "ordinality")

d %>% 
  filter(year > 0) %>%
  gather(measure, value, deviance, linr2, ans, ordinality) %>%
  group_by(year, measure) %>% 
  do(tidy(t.test(.$value[.$condition == "abacus"], 
                 .$value[.$condition == "control"])))
## Source: local data frame [12 x 10]
## Groups: year, measure
## 
##    year    measure     estimate estimate1 estimate2  statistic    p.value
## 1     1   deviance -0.080161367 0.6653488 0.7455102 -1.6349326 0.10385750
## 2     1      linr2  0.015889891 0.3747674 0.3588776  0.5947372 0.55279248
## 3     1        ans -0.002817625 0.1758140 0.1786316 -0.1979361 0.84333389
## 4     1 ordinality  0.021198386 0.8089535 0.7877551  2.1472626 0.03312995
## 5     2   deviance -0.117626263 0.6504545 0.7680808 -2.2387863 0.02637566
## 6     2      linr2  0.010303030 0.3572727 0.3469697  0.3780748 0.70582175
## 7     2        ans -0.002969775 0.1446591 0.1476289 -0.3397382 0.73446692
## 8     2 ordinality  0.014166667 0.8011364 0.7869697  1.3926937 0.16539755
## 9     3   deviance -0.072769070 0.6596552 0.7324242 -1.4796581 0.14067668
## 10    3      linr2  0.067293626 0.3918391 0.3245455  2.3208574 0.02146247
## 11    3        ans -0.011578700 0.1322989 0.1438776 -1.1261363 0.26159050
## 12    3 ordinality  0.011323581 0.8042529 0.7929293  1.1191883 0.26453958
## Variables not shown: parameter (dbl), conf.low (dbl), conf.high (dbl)

PLOTS

ANS - Without intervention split

qplot(ans, standardized, facets=~year, 
      data=d) + 
  geom_smooth(method="lm") +
  ylab("Standardized Test Composite") +
  xlab("Weber Fraction")        

And with:

qplot(ans, standardized, facets=~year, 
      col=condition,
      data=d) + 
  geom_smooth(method="lm") +
  ylab("Standardized Test Composite") +
  xlab("Weber Fraction")        

Now do the same for deviance: ```{/r, warnings=FALSE} qplot(deviance, standardized, facets=~year, data=d) + geom_smooth(method=“lm”) + ylab(“Standardized Test Composite”) + xlab(“PAE”)

qplot(deviance, standardized, facets=~year, col=condition, data=d) + geom_smooth(method=“lm”) + ylab(“Standardized Test Composite”) + xlab(“PAE”)
```

and for linear \(r^2\):

qplot(linr2, standardized, facets=~year,
      data=d) + 
  geom_smooth(method="lm") +
  ylab("Standardized Test Composite") +
  xlab("Linear r^2")    

qplot(linr2, standardized, facets=~year, 
      col=condition,
      data=d) + 
  geom_smooth(method="lm") +
  ylab("Standardized Test Composite") +
  xlab("Linear r^2")       

and for ordinality:

qplot(ordinality, standardized, facets=~year,
      data=d) + 
  geom_smooth(method="lm") +
  ylab("Standardized Test Composite") +
  xlab("Ordinality")    

qplot(ordinality, standardized, facets=~year, 
      col=condition,
      data=d) + 
  geom_smooth(method="lm") +
  ylab("Standardized Test Composite") +
  xlab("Ordinality")       

Models

Standardize predictors.

# this is insane, but PCA will choke if the scaled variables have included in them the scaling attributes (and this is what R does by default). 
sscale <- function (x) {as.numeric(scale(x))}

md <- d %>% 
  gather(measure, value, deviance, linr2, ans, ordinality) %>%  
  group_by(year, measure) %>%
  mutate(standardized.scale = sscale(standardized),
         value.scale = sscale(value),
         wiat.scale = sscale(wiat), 
         woodcock.scale = sscale(woodcock),
         math.scale = sscale(math), 
         placeval.scale = sscale(placeval),
         arith.scale = sscale(arith), 
         mental.rot.scale = sscale(mental.rot),
         verbalwm.scale = sscale(verbalwm),
         spatialwm.scale = sscale(spatialwm),
         ravens.scale = sscale(ravens)) 

Now fit models. This uses the broom package - really very nifty.

NOTE, these are scaled values predicting scaled values.

std.baseline.models <- md %>%
  group_by(year, measure) %>%
  filter(!(year == 0 & measure != "ans")) %>%
  do(tidy(lm(standardized.scale ~ value.scale, data = .))) %>%
  filter(term != "(Intercept)") %>%
  data.frame
std.baseline.models
##    year    measure        term      estimate  std.error   statistic
## 1     0        ans value.scale -0.2861136935 0.07585437 -3.77188164
## 2     1   deviance value.scale -0.0736166477 0.07429885 -0.99081808
## 3     1      linr2 value.scale  0.1031406510 0.07410564  1.39180566
## 4     1        ans value.scale -0.1210668281 0.07315981 -1.65482699
## 5     1 ordinality value.scale  0.0008627001 0.07449894  0.01158003
## 6     2   deviance value.scale -0.1447405155 0.07301173 -1.98242838
## 7     2      linr2 value.scale  0.1910128279 0.07337888  2.60310359
## 8     2        ans value.scale -0.1049879875 0.07382586 -1.42210317
## 9     2 ordinality value.scale  0.2081562771 0.07257958  2.86797306
## 10    3   deviance value.scale -0.2283394329 0.07187531 -3.17688289
## 11    3      linr2 value.scale  0.2773777821 0.07093151  3.91050137
## 12    3        ans value.scale -0.2259644860 0.07165180 -3.15364723
## 13    3 ordinality value.scale  0.2275793284 0.07188841  3.16573055
##         p.value
## 1  0.0002276970
## 2  0.3230898640
## 3  0.1656798389
## 4  0.0997119177
## 5  0.9907733616
## 6  0.0489197083
## 7  0.0099918231
## 8  0.1567077999
## 9  0.0046136812
## 10 0.0017460976
## 11 0.0001294263
## 12 0.0018849401
## 13 0.0018107944

Check year 3 correlations

cor.test(d$standardized[d$year==3],d$linr2[d$year==3])
## 
##  Pearson's product-moment correlation
## 
## data:  d$standardized[d$year == 3] and d$linr2[d$year == 3]
## t = 3.9105, df = 184, p-value = 0.0001294
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1386516 0.4047527
## sample estimates:
##       cor 
## 0.2770049
cor.test(d$standardized[d$year==3],d$ans[d$year==3])
## 
##  Pearson's product-moment correlation
## 
## data:  d$standardized[d$year == 3] and d$ans[d$year == 3]
## t = -3.1536, df = 183, p-value = 0.001885
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.35952906 -0.08557076
## sample estimates:
##        cor 
## -0.2270366
cor.test(d$standardized[d$year==3],d$deviance[d$year==3])
## 
##  Pearson's product-moment correlation
## 
## data:  d$standardized[d$year == 3] and d$deviance[d$year == 3]
## t = -3.1769, df = 184, p-value = 0.001746
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.36009736 -0.08700778
## sample estimates:
##        cor 
## -0.2280325
cor.test(d$standardized[d$year==3],d$ordinality[d$year==3])
## 
##  Pearson's product-moment correlation
## 
## data:  d$standardized[d$year == 3] and d$ordinality[d$year == 3]
## t = 3.1657, df = 184, p-value = 0.001811
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.08621322 0.35940039
## sample estimates:
##       cor 
## 0.2272734

And here are the full models.

std.models <- md %>%
  group_by(year, measure) %>%
  filter(!(year == 0 & measure != "ans")) %>%
  do(tidy(lm(standardized.scale ~ value.scale + 
               mental.rot.scale + 
               verbalwm.scale + 
               spatialwm.scale + 
               ravens.scale + 
               age + condition, data = .)))

options(dplyr.width = Inf)
print.data.frame(std.models)
##     year    measure             term      estimate  std.error    statistic
## 1      0        ans      (Intercept) -0.7856538941 0.97930903 -0.802253295
## 2      0        ans      value.scale -0.2382604904 0.07869646 -3.027588537
## 3      0        ans mental.rot.scale -0.1506849858 0.08279761 -1.819919451
## 4      0        ans   verbalwm.scale  0.1827323959 0.09210359  1.983987871
## 5      0        ans  spatialwm.scale  0.2298199171 0.09478945  2.424530612
## 6      0        ans     ravens.scale -0.0032556306 0.08233557 -0.039540997
## 7      0        ans              age  0.1262516923 0.14629239  0.863009287
## 8      0        ans  conditionabacus  0.0305031720 0.15180823  0.200932275
## 9      1   deviance      (Intercept)  0.0704332170 0.87808443  0.080212352
## 10     1   deviance      value.scale  0.0069503802 0.07071760  0.098283598
## 11     1   deviance mental.rot.scale  0.1684011132 0.08202233  2.053112944
## 12     1   deviance   verbalwm.scale  0.1168479389 0.07678948  1.521665910
## 13     1   deviance  spatialwm.scale  0.0875094913 0.07672455  1.140566990
## 14     1   deviance     ravens.scale  0.2263509344 0.08458085  2.676148647
## 15     1   deviance              age -0.0304643051 0.13145994 -0.231738320
## 16     1   deviance  conditionabacus  0.2868962144 0.14741978  1.946117463
## 17     1      linr2      (Intercept)  0.0660351616 0.87843933  0.075173275
## 18     1      linr2      value.scale  0.0089746826 0.07048610  0.127325562
## 19     1      linr2 mental.rot.scale  0.1669835035 0.08228118  2.029425182
## 20     1      linr2   verbalwm.scale  0.1149181433 0.07730003  1.486650723
## 21     1      linr2  spatialwm.scale  0.0858767514 0.07573664  1.133886506
## 22     1      linr2     ravens.scale  0.2274777970 0.08449940  2.692064207
## 23     1      linr2              age -0.0297437056 0.13150432 -0.226180453
## 24     1      linr2  conditionabacus  0.2862859782 0.14722525  1.944544074
## 25     1        ans      (Intercept) -0.1421182708 0.92169939 -0.154191564
## 26     1        ans      value.scale -0.0848020275 0.06879832 -1.232617643
## 27     1        ans mental.rot.scale  0.1392803795 0.08123220  1.714595685
## 28     1        ans   verbalwm.scale  0.1313498236 0.07610603  1.725879371
## 29     1        ans  spatialwm.scale  0.0851861084 0.07508451  1.134536400
## 30     1        ans     ravens.scale  0.2132227330 0.08329833  2.559748108
## 31     1        ans              age  0.0095505003 0.13858947  0.068912163
## 32     1        ans  conditionabacus  0.2263160447 0.14656650  1.544118489
## 33     1 ordinality      (Intercept)  0.0876727683 0.87616117  0.100064659
## 34     1 ordinality      value.scale -0.0646030826 0.07164768 -0.901677256
## 35     1 ordinality mental.rot.scale  0.1689033636 0.08172022  2.066849190
## 36     1 ordinality   verbalwm.scale  0.1263959819 0.07720328  1.637184133
## 37     1 ordinality  spatialwm.scale  0.0991810438 0.07684121  1.290727203
## 38     1 ordinality     ravens.scale  0.2179444882 0.08476923  2.571033043
## 39     1 ordinality              age -0.0335287511 0.13117888 -0.255595660
## 40     1 ordinality  conditionabacus  0.2931081377 0.14707403  1.992929304
## 41     2   deviance      (Intercept) -0.5018657569 0.88948097 -0.564223155
## 42     2   deviance      value.scale -0.0511800660 0.07373884 -0.694071980
## 43     2   deviance mental.rot.scale  0.1230791522 0.08065492  1.525996802
## 44     2   deviance   verbalwm.scale  0.1281293547 0.07541566  1.698975429
## 45     2   deviance  spatialwm.scale -0.0340836956 0.07332083 -0.464856916
## 46     2   deviance     ravens.scale  0.2693148811 0.08271429  3.255965458
## 47     2   deviance              age  0.0580492047 0.13358344  0.434553904
## 48     2   deviance  conditionabacus  0.2639209643 0.14321187  1.842870793
## 49     2      linr2      (Intercept) -0.3876979388 0.87850525 -0.441315448
## 50     2      linr2      value.scale  0.1380805185 0.07299928  1.891532724
## 51     2      linr2 mental.rot.scale  0.1124738978 0.08013608  1.403536268
## 52     2      linr2   verbalwm.scale  0.1099306123 0.07521069  1.461635399
## 53     2      linr2  spatialwm.scale -0.0389742974 0.07268825 -0.536184256
## 54     2      linr2     ravens.scale  0.2806527199 0.08098731  3.465391137
## 55     2      linr2              age  0.0396760723 0.13179446  0.301045078
## 56     2      linr2  conditionabacus  0.2764312569 0.14042647  1.968512480
## 57     2        ans      (Intercept) -0.3072567066 0.94018041 -0.326806113
## 58     2        ans      value.scale -0.0301034269 0.07696909 -0.391110614
## 59     2        ans mental.rot.scale  0.1146072990 0.08415605  1.361842738
## 60     2        ans   verbalwm.scale  0.1365010632 0.07658785  1.782280867
## 61     2        ans  spatialwm.scale -0.0410087042 0.07553182 -0.542932823
## 62     2        ans     ravens.scale  0.2833881112 0.08283879  3.420958934
## 63     2        ans              age  0.0270128435 0.14179116  0.190511472
## 64     2        ans  conditionabacus  0.2836986921 0.14350583  1.976914100
## 65     2 ordinality      (Intercept) -0.7205916314 0.88609373 -0.813222806
## 66     2 ordinality      value.scale  0.1512482795 0.07301202  2.071553260
## 67     2 ordinality mental.rot.scale  0.1140304473 0.07989318  1.427286430
## 68     2 ordinality   verbalwm.scale  0.1137262272 0.07457634  1.524963982
## 69     2 ordinality  spatialwm.scale -0.0273514290 0.07247060 -0.377414128
## 70     2 ordinality     ravens.scale  0.2668443782 0.08099826  3.294445631
## 71     2 ordinality              age  0.0925382873 0.13313316  0.695080665
## 72     2 ordinality  conditionabacus  0.2464247627 0.14098063  1.747933438
## 73     3   deviance      (Intercept) -0.2560927282 0.80749465 -0.317144799
## 74     3   deviance      value.scale -0.0357769329 0.06927171 -0.516472504
## 75     3   deviance mental.rot.scale  0.2709432040 0.07600035  3.565025514
## 76     3   deviance   verbalwm.scale  0.1099748891 0.06626194  1.659699161
## 77     3   deviance  spatialwm.scale  0.0994894786 0.07112983  1.398702634
## 78     3   deviance     ravens.scale  0.2433194324 0.07590442  3.205603074
## 79     3   deviance              age  0.0256576208 0.12083417  0.212337463
## 80     3   deviance  conditionabacus  0.1907479142 0.13077011  1.458650670
## 81     3      linr2      (Intercept) -0.0655470230 0.80968285 -0.080953948
## 82     3      linr2      value.scale  0.1128325528 0.06946138  1.624392521
## 83     3      linr2 mental.rot.scale  0.2427626639 0.07753767  3.130899542
## 84     3      linr2   verbalwm.scale  0.0959844742 0.06643793  1.444724029
## 85     3      linr2  spatialwm.scale  0.1108330877 0.06911768  1.603541794
## 86     3      linr2     ravens.scale  0.2449441035 0.07497269  3.267110975
## 87     3      linr2              age -0.0008532564 0.12090824 -0.007057057
## 88     3      linr2  conditionabacus  0.1614555244 0.13088925  1.233527735
## 89     3        ans      (Intercept) -0.3544991433 0.83241480 -0.425868383
## 90     3        ans      value.scale -0.0629947326 0.06825250 -0.922965986
## 91     3        ans mental.rot.scale  0.2684758022 0.07572026  3.545627172
## 92     3        ans   verbalwm.scale  0.1039107589 0.06691078  1.552974742
## 93     3        ans  spatialwm.scale  0.0879826488 0.07247507  1.213971244
## 94     3        ans     ravens.scale  0.2384733210 0.07604512  3.135944900
## 95     3        ans              age  0.0412869395 0.12490245  0.330553476
## 96     3        ans  conditionabacus  0.1844731761 0.13059066  1.412606203
## 97     3 ordinality      (Intercept) -0.2044912280 0.80854880 -0.252911424
## 98     3 ordinality      value.scale  0.0528422831 0.06870094  0.769163866
## 99     3 ordinality mental.rot.scale  0.2691781146 0.07574292  3.553838830
## 100    3 ordinality   verbalwm.scale  0.1073684321 0.06635353  1.618126885
## 101    3 ordinality  spatialwm.scale  0.0977280214 0.07056627  1.384911340
## 102    3 ordinality     ravens.scale  0.2416496514 0.07578199  3.188747579
## 103    3 ordinality              age  0.0180278107 0.12089009  0.149125624
## 104    3 ordinality  conditionabacus  0.1896118295 0.13017170  1.456628716
##          p.value
## 1   0.4237958964
## 2   0.0029461438
## 3   0.0709542685
## 4   0.0492556845
## 5   0.0166316159
## 6   0.9685166315
## 7   0.3896406260
## 8   0.8410495331
## 9   0.9361644093
## 10  0.9218249689
## 11  0.0416214241
## 12  0.1299840292
## 13  0.2556832131
## 14  0.0081894521
## 15  0.8170249786
## 16  0.0533190830
## 17  0.9401668218
## 18  0.8988359129
## 19  0.0440026059
## 20  0.1389926135
## 21  0.2584672279
## 22  0.0078244039
## 23  0.8213376636
## 24  0.0535099277
## 25  0.8776482769
## 26  0.2194825325
## 27  0.0883083896
## 28  0.0862525322
## 29  0.2582251573
## 30  0.0113779623
## 31  0.9451434874
## 32  0.1244870632
## 33  0.9204130166
## 34  0.3685269175
## 35  0.0402916937
## 36  0.1034749939
## 37  0.1985823287
## 38  0.0110117305
## 39  0.7985776068
## 40  0.0478980056
## 41  0.5733459095
## 42  0.4885844316
## 43  0.1288689680
## 44  0.0911530623
## 45  0.6426287760
## 46  0.0013640044
## 47  0.6644373417
## 48  0.0670893756
## 49  0.6595454225
## 50  0.0602538121
## 51  0.1622809782
## 52  0.1456871485
## 53  0.5925321054
## 54  0.0006702930
## 55  0.7637480085
## 56  0.0506350122
## 57  0.7442213401
## 58  0.6962110432
## 59  0.1750708691
## 60  0.0765095661
## 61  0.5878956429
## 62  0.0007835058
## 63  0.8491383711
## 64  0.0496877453
## 65  0.4172278081
## 66  0.0398175172
## 67  0.1553315807
## 68  0.1291262918
## 69  0.7063366207
## 70  0.0012000617
## 71  0.4879537205
## 72  0.0822809528
## 73  0.7515161450
## 74  0.6061839930
## 75  0.0004705116
## 76  0.0987866058
## 77  0.1636923296
## 78  0.0016053453
## 79  0.8320937454
## 80  0.1464744448
## 81  0.9355721535
## 82  0.1061126096
## 83  0.0020463495
## 84  0.1503441964
## 85  0.1106392139
## 86  0.0013103222
## 87  0.9943774607
## 88  0.2190522339
## 89  0.6707361801
## 90  0.3573178456
## 91  0.0005046884
## 92  0.1222672477
## 93  0.2264229451
## 94  0.0020151761
## 95  0.7413840986
## 96  0.1595784209
## 97  0.8006369287
## 98  0.4428448800
## 99  0.0004895490
## 100 0.1074570113
## 101 0.1678629307
## 102 0.0016963512
## 103 0.8816282024
## 104 0.1470314608
gelb<-as.numeric(std.models$p.value)
small<-subset(std.models, gelb<.05)
print.data.frame(small)
##    year    measure             term   estimate  std.error statistic
## 1     0        ans      value.scale -0.2382605 0.07869646 -3.027589
## 2     0        ans   verbalwm.scale  0.1827324 0.09210359  1.983988
## 3     0        ans  spatialwm.scale  0.2298199 0.09478945  2.424531
## 4     1   deviance mental.rot.scale  0.1684011 0.08202233  2.053113
## 5     1   deviance     ravens.scale  0.2263509 0.08458085  2.676149
## 6     1      linr2 mental.rot.scale  0.1669835 0.08228118  2.029425
## 7     1      linr2     ravens.scale  0.2274778 0.08449940  2.692064
## 8     1        ans     ravens.scale  0.2132227 0.08329833  2.559748
## 9     1 ordinality mental.rot.scale  0.1689034 0.08172022  2.066849
## 10    1 ordinality     ravens.scale  0.2179445 0.08476923  2.571033
## 11    1 ordinality  conditionabacus  0.2931081 0.14707403  1.992929
## 12    2   deviance     ravens.scale  0.2693149 0.08271429  3.255965
## 13    2      linr2     ravens.scale  0.2806527 0.08098731  3.465391
## 14    2        ans     ravens.scale  0.2833881 0.08283879  3.420959
## 15    2        ans  conditionabacus  0.2836987 0.14350583  1.976914
## 16    2 ordinality      value.scale  0.1512483 0.07301202  2.071553
## 17    2 ordinality     ravens.scale  0.2668444 0.08099826  3.294446
## 18    3   deviance mental.rot.scale  0.2709432 0.07600035  3.565026
## 19    3   deviance     ravens.scale  0.2433194 0.07590442  3.205603
## 20    3      linr2 mental.rot.scale  0.2427627 0.07753767  3.130900
## 21    3      linr2     ravens.scale  0.2449441 0.07497269  3.267111
## 22    3        ans mental.rot.scale  0.2684758 0.07572026  3.545627
## 23    3        ans     ravens.scale  0.2384733 0.07604512  3.135945
## 24    3 ordinality mental.rot.scale  0.2691781 0.07574292  3.553839
## 25    3 ordinality     ravens.scale  0.2416497 0.07578199  3.188748
##         p.value
## 1  0.0029461438
## 2  0.0492556845
## 3  0.0166316159
## 4  0.0416214241
## 5  0.0081894521
## 6  0.0440026059
## 7  0.0078244039
## 8  0.0113779623
## 9  0.0402916937
## 10 0.0110117305
## 11 0.0478980056
## 12 0.0013640044
## 13 0.0006702930
## 14 0.0007835058
## 15 0.0496877453
## 16 0.0398175172
## 17 0.0012000617
## 18 0.0004705116
## 19 0.0016053453
## 20 0.0020463495
## 21 0.0013103222
## 22 0.0005046884
## 23 0.0020151761
## 24 0.0004895490
## 25 0.0016963512

And make pretty plot.

std.models$term <- factor(std.models$term, 
                           levels = c("(Intercept)",
                                      "value.scale", "mental.rot.scale",
                                      "spatialwm.scale", "verbalwm.scale",
                                      "ravens.scale","age", "conditionabacus"), 
                           labels = c("Intercept", 
                                      "Predictor", "Mental Rotation",
                                      "Spatial WM", "Verbal WM",
                                      "Raven's", "Age", "Intervention"))

std.models$measure <- factor(std.models$measure,
                              levels = c("ans","deviance","linr2","ordinality"),
                              labels = c("ANS","PAE","Linear r^2", 
                                         "Ordinality"))
  
qplot(term, estimate, 
      fill = term, 
      ymin = estimate - std.error, ymax = estimate + std.error,
      geom = c("bar", "linerange"), stat = "identity", 
      facets = measure ~ year, 
      data=filter(std.models, term != "Intercept")) + 
  geom_hline(yintercept=0, lty=2) + 
  xlab("Predictor") + 
  ylab("Standardized Beta Weight") + 
  scale_fill_discrete(name="Predictor") +
  theme(axis.text.x = element_text(angle = 90, vjust=.5, hjust = 1)) 

Add PCA to the dataset. Note that PC1 is flipped just to make life easier.

# need this to get around standard evals
pc1 <- function(x,y,z,m,n) {
  sscale(as.numeric(prcomp(~x + y + z + m + n)$x[,1]))
}

# need to filter to complete cases
pmd <- md %>%
  filter(complete.cases(wiat.scale, woodcock.scale, math.scale, 
                        placeval.scale, arith.scale)) %>%
  mutate(pc1 = -pc1(wiat.scale, woodcock.scale, math.scale, 
                   placeval.scale, arith.scale))

PCA Models and plot.

pca.baseline.models <- pmd %>%
  group_by(year, measure) %>%
  filter(!(year == 0 & measure != "ans")) %>%
  do(tidy(lm(pc1 ~ value.scale, data = .))) %>%
  filter(term != "(Intercept)") %>%
  data.frame
pca.baseline.models
##    year    measure        term   estimate  std.error statistic
## 1     0        ans value.scale  0.2960511 0.07929057  3.733749
## 2     1   deviance value.scale -0.1676604 0.07346686 -2.282123
## 3     1      linr2 value.scale  0.1581167 0.07354619  2.149897
## 4     1        ans value.scale -0.1294271 0.07288281 -1.775825
## 5     1 ordinality value.scale  0.1054681 0.07430609  1.419374
## 6     2   deviance value.scale -0.1873279 0.07311760 -2.562009
## 7     2      linr2 value.scale  0.2497652 0.07302858  3.420102
## 8     2        ans value.scale -0.1474893 0.07411458 -1.990018
## 9     2 ordinality value.scale  0.1914047 0.07492045  2.554772
## 10    3   deviance value.scale -0.2476530 0.07146468 -3.465390
## 11    3      linr2 value.scale  0.3378703 0.06971814  4.846232
## 12    3        ans value.scale -0.2617409 0.07125388 -3.673356
## 13    3 ordinality value.scale  0.2735804 0.07165864  3.817828
##         p.value
## 1  2.794896e-04
## 2  2.364572e-02
## 3  3.289062e-02
## 4  7.747032e-02
## 5  1.575100e-01
## 6  1.122584e-02
## 7  7.744574e-04
## 8  4.811134e-02
## 9  1.145338e-02
## 10 6.603674e-04
## 11 2.684857e-06
## 12 3.150773e-04
## 13 1.844388e-04
pca.models <- pmd %>%
  group_by(year, measure) %>%
  filter(!(year == 0 & measure != "ans")) %>%
  do(tidy(lm(pc1 ~ value.scale + 
               mental.rot.scale + 
               verbalwm.scale + 
               spatialwm.scale + 
               ravens.scale + 
               age + condition, data = .)))

options(dplyr.width = Inf)

print.data.frame(pca.models)
##     year    measure             term     estimate  std.error   statistic
## 1      0        ans      (Intercept)  1.497922211 0.97045769  1.54352140
## 2      0        ans      value.scale  0.185976822 0.07924519  2.34685317
## 3      0        ans mental.rot.scale  0.072935906 0.08421307  0.86608773
## 4      0        ans   verbalwm.scale -0.202937260 0.08948352 -2.26787308
## 5      0        ans  spatialwm.scale -0.261960198 0.09449228 -2.77229205
## 6      0        ans     ravens.scale -0.114758661 0.08801196 -1.30389846
## 7      0        ans              age -0.229856115 0.14462246 -1.58935285
## 8      0        ans  conditionabacus -0.093480549 0.15332473 -0.60968995
## 9      1   deviance      (Intercept) -0.393256073 0.86642719 -0.45388243
## 10     1   deviance      value.scale -0.076369306 0.06956242 -1.09785292
## 11     1   deviance mental.rot.scale  0.137902965 0.08062877  1.71034434
## 12     1   deviance   verbalwm.scale  0.154817774 0.07547741  2.05118013
## 13     1   deviance  spatialwm.scale  0.137621074 0.07537677  1.82577580
## 14     1   deviance     ravens.scale  0.215687795 0.08310873  2.59524824
## 15     1   deviance              age  0.041116779 0.12989342  0.31654243
## 16     1   deviance  conditionabacus  0.260282871 0.14542286  1.78983458
## 17     1      linr2      (Intercept) -0.411056655 0.86683764 -0.47420259
## 18     1      linr2      value.scale  0.076027446 0.06928729  1.09727842
## 19     1      linr2 mental.rot.scale  0.134040458 0.08088218  1.65723111
## 20     1      linr2   verbalwm.scale  0.149613179 0.07597723  1.96918457
## 21     1      linr2  spatialwm.scale  0.147930968 0.07442020  1.98777981
## 22     1      linr2     ravens.scale  0.214289364 0.08303426  2.58073440
## 23     1      linr2              age  0.043047989 0.12994334  0.33128276
## 24     1      linr2  conditionabacus  0.269851367 0.14519564  1.85853635
## 25     1        ans      (Intercept) -0.719905636 0.90360982 -0.79669966
## 26     1        ans      value.scale -0.091798116 0.06736004 -1.36279788
## 27     1        ans mental.rot.scale  0.110370490 0.07934901  1.39094976
## 28     1        ans   verbalwm.scale  0.179227511 0.07434468  2.41076429
## 29     1        ans  spatialwm.scale  0.148767080 0.07333249  2.02866532
## 30     1        ans     ravens.scale  0.193001655 0.08136940  2.37191943
## 31     1        ans              age  0.099131013 0.13604215  0.72867867
## 32     1        ans  conditionabacus  0.196803424 0.14371770  1.36937501
## 33     1 ordinality      (Intercept) -0.387559399 0.86962288 -0.44566375
## 34     1 ordinality      value.scale  0.023544324 0.07111527  0.33107270
## 35     1 ordinality mental.rot.scale  0.142111839 0.08079549  1.75890797
## 36     1 ordinality   verbalwm.scale  0.157371241 0.07629260  2.06273283
## 37     1 ordinality  spatialwm.scale  0.146513481 0.07591240  1.93003354
## 38     1 ordinality     ravens.scale  0.212816809 0.08375460  2.54095678
## 39     1 ordinality              age  0.039723398 0.13038597  0.30466006
## 40     1 ordinality  conditionabacus  0.266280511 0.14594540  1.82452147
## 41     2   deviance      (Intercept) -0.559636755 0.87668416 -0.63835618
## 42     2   deviance      value.scale -0.072862999 0.07099878 -1.02625711
## 43     2   deviance mental.rot.scale  0.105192762 0.07799517  1.34870861
## 44     2   deviance   verbalwm.scale  0.171976503 0.07220109  2.38191012
## 45     2   deviance  spatialwm.scale -0.012544974 0.07086466 -0.17702723
## 46     2   deviance     ravens.scale  0.301492650 0.07928570  3.80261067
## 47     2   deviance              age  0.055304097 0.13162517  0.42016355
## 48     2   deviance  conditionabacus  0.414240079 0.13779835  3.00613242
## 49     2      linr2      (Intercept) -0.440782477 0.85431608 -0.51594777
## 50     2      linr2      value.scale  0.191623271 0.06961279  2.75270222
## 51     2      linr2 mental.rot.scale  0.089100273 0.07673112  1.16120132
## 52     2      linr2   verbalwm.scale  0.145760406 0.07140294  2.04137808
## 53     2      linr2  spatialwm.scale -0.018563201 0.06949943 -0.26709862
## 54     2      linr2     ravens.scale  0.318710012 0.07676868  4.15156282
## 55     2      linr2              age  0.035902866 0.12811956  0.28022940
## 56     2      linr2  conditionabacus  0.433918051 0.13400545  3.23806259
## 57     2        ans      (Intercept) -0.311710269 0.91191696 -0.34181870
## 58     2        ans      value.scale -0.071477985 0.07407185 -0.96498183
## 59     2        ans mental.rot.scale  0.086794558 0.08049890  1.07820795
## 60     2        ans   verbalwm.scale  0.174039099 0.07312008  2.38018203
## 61     2        ans  spatialwm.scale -0.020299232 0.07245154 -0.28017670
## 62     2        ans     ravens.scale  0.318265539 0.07901458  4.02793447
## 63     2        ans              age  0.016343796 0.13727361  0.11906000
## 64     2        ans  conditionabacus  0.434433459 0.13760369  3.15713528
## 65     2 ordinality      (Intercept) -0.691272223 0.87443807 -0.79053308
## 66     2 ordinality      value.scale  0.130566165 0.07138147  1.82913258
## 67     2 ordinality mental.rot.scale  0.101307039 0.07749056  1.30734685
## 68     2 ordinality   verbalwm.scale  0.164190252 0.07171437  2.28950280
## 69     2 ordinality  spatialwm.scale -0.003238722 0.07031302 -0.04606147
## 70     2 ordinality     ravens.scale  0.305298433 0.07787643  3.92029326
## 71     2 ordinality              age  0.075224262 0.13125771  0.57310357
## 72     2 ordinality  conditionabacus  0.411803017 0.13615714  3.02446873
## 73     3   deviance      (Intercept)  0.201030221 0.77841509  0.25825581
## 74     3   deviance      value.scale -0.033882164 0.06661796 -0.50860401
## 75     3   deviance mental.rot.scale  0.265338664 0.07480796  3.54693091
## 76     3   deviance   verbalwm.scale  0.093297204 0.06367089  1.46530389
## 77     3   deviance  spatialwm.scale  0.133617372 0.06872634  1.94419458
## 78     3   deviance     ravens.scale  0.288460659 0.07308947  3.94667858
## 79     3   deviance              age -0.053056682 0.11630778 -0.45617481
## 80     3   deviance  conditionabacus  0.328513763 0.12625532  2.60197950
## 81     3      linr2      (Intercept)  0.488534014 0.77079205  0.63380780
## 82     3      linr2      value.scale  0.174439128 0.06581358  2.65050345
## 83     3      linr2 mental.rot.scale  0.219014358 0.07509538  2.91648235
## 84     3      linr2   verbalwm.scale  0.070574026 0.06304830  1.11936438
## 85     3      linr2  spatialwm.scale  0.146887808 0.06588906  2.22932000
## 86     3      linr2     ravens.scale  0.288150617 0.07133822  4.03921794
## 87     3      linr2              age -0.092625073 0.11493704 -0.80587665
## 88     3      linr2  conditionabacus  0.278366611 0.12466335  2.23294668
## 89     3        ans      (Intercept)  0.112100686 0.79960250  0.14019552
## 90     3        ans      value.scale -0.091779712 0.06539806 -1.40340108
## 91     3        ans mental.rot.scale  0.260719124 0.07407022  3.51989114
## 92     3        ans   verbalwm.scale  0.086171485 0.06410000  1.34432891
## 93     3        ans  spatialwm.scale  0.118246779 0.06969060  1.69673923
## 94     3        ans     ravens.scale  0.279641849 0.07308300  3.82635973
## 95     3        ans              age -0.038693255 0.11982097 -0.32292558
## 96     3        ans  conditionabacus  0.318335301 0.12554495  2.53562813
## 97     3 ordinality      (Intercept)  0.283253764 0.77575484  0.36513309
## 98     3 ordinality      value.scale  0.097249409 0.06583697  1.47712470
## 99     3 ordinality mental.rot.scale  0.258919954 0.07393976  3.50176886
## 100    3 ordinality   verbalwm.scale  0.086623769 0.06349026  1.36436311
## 101    3 ordinality  spatialwm.scale  0.124590058 0.06769418  1.84048404
## 102    3 ordinality     ravens.scale  0.280743072 0.07272759  3.86020023
## 103    3 ordinality              age -0.064727692 0.11582482 -0.55884127
## 104    3 ordinality  conditionabacus  0.319673573 0.12493894  2.55863851
##          p.value
## 1   1.254758e-01
## 2   2.065866e-02
## 3   3.882612e-01
## 4   2.522071e-02
## 5   6.503401e-03
## 6   1.948950e-01
## 7   1.147503e-01
## 8   5.432805e-01
## 9   6.505063e-01
## 10  2.738586e-01
## 11  8.907046e-02
## 12  4.182095e-02
## 13  6.968097e-02
## 14  1.029879e-02
## 15  7.519885e-01
## 16  7.530337e-02
## 17  6.359786e-01
## 18  2.741089e-01
## 19  9.936150e-02
## 20  5.059636e-02
## 21  4.848038e-02
## 22  1.072408e-02
## 23  7.408487e-01
## 24  6.486347e-02
## 25  4.267844e-01
## 26  1.748260e-01
## 27  1.661365e-01
## 28  1.703141e-02
## 29  4.411998e-02
## 30  1.886275e-02
## 31  4.672439e-01
## 32  1.727658e-01
## 33  6.564208e-01
## 34  7.410071e-01
## 35  8.043567e-02
## 36  4.069567e-02
## 37  5.530750e-02
## 38  1.197156e-02
## 39  7.610067e-01
## 40  6.987117e-02
## 41  5.241215e-01
## 42  3.062635e-01
## 43  1.792677e-01
## 44  1.835480e-02
## 45  8.597028e-01
## 46  2.008624e-04
## 47  6.749093e-01
## 48  3.057078e-03
## 49  6.065774e-01
## 50  6.568543e-03
## 51  2.472271e-01
## 52  4.279644e-02
## 53  7.897248e-01
## 54  5.267619e-05
## 55  7.796504e-01
## 56  1.453375e-03
## 57  7.329223e-01
## 58  3.359661e-01
## 59  2.825147e-01
## 60  1.844511e-02
## 61  7.796929e-01
## 62  8.567010e-05
## 63  9.053727e-01
## 64  1.894650e-03
## 65  4.303444e-01
## 66  6.917410e-02
## 67  1.929028e-01
## 68  2.330961e-02
## 69  9.633166e-01
## 70  1.291290e-04
## 71  5.673501e-01
## 72  2.886892e-03
## 73  7.965201e-01
## 74  6.116854e-01
## 75  5.030810e-04
## 76  1.446743e-01
## 77  5.351301e-02
## 78  1.155482e-04
## 79  6.488431e-01
## 80  1.008181e-02
## 81  5.270531e-01
## 82  8.791686e-03
## 83  4.015162e-03
## 84  2.645546e-01
## 85  2.709501e-02
## 86  8.086261e-05
## 87  4.214332e-01
## 88  2.685006e-02
## 89  8.886715e-01
## 90  1.623212e-01
## 91  5.541413e-04
## 92  1.806332e-01
## 93  9.157615e-02
## 94  1.824602e-04
## 95  7.471485e-01
## 96  1.212612e-02
## 97  7.154633e-01
## 98  1.414815e-01
## 99  5.896754e-04
## 100 1.742462e-01
## 101 6.743037e-02
## 102 1.604113e-04
## 103 5.770011e-01
## 104 1.137530e-02
gelb2<-as.numeric(pca.models$p.value)
small2<-subset(pca.models, gelb2<.05)
print.data.frame(small2)
##    year    measure             term   estimate  std.error statistic
## 1     0        ans      value.scale  0.1859768 0.07924519  2.346853
## 2     0        ans   verbalwm.scale -0.2029373 0.08948352 -2.267873
## 3     0        ans  spatialwm.scale -0.2619602 0.09449228 -2.772292
## 4     1   deviance   verbalwm.scale  0.1548178 0.07547741  2.051180
## 5     1   deviance     ravens.scale  0.2156878 0.08310873  2.595248
## 6     1      linr2  spatialwm.scale  0.1479310 0.07442020  1.987780
## 7     1      linr2     ravens.scale  0.2142894 0.08303426  2.580734
## 8     1        ans   verbalwm.scale  0.1792275 0.07434468  2.410764
## 9     1        ans  spatialwm.scale  0.1487671 0.07333249  2.028665
## 10    1        ans     ravens.scale  0.1930017 0.08136940  2.371919
## 11    1 ordinality   verbalwm.scale  0.1573712 0.07629260  2.062733
## 12    1 ordinality     ravens.scale  0.2128168 0.08375460  2.540957
## 13    2   deviance   verbalwm.scale  0.1719765 0.07220109  2.381910
## 14    2   deviance     ravens.scale  0.3014926 0.07928570  3.802611
## 15    2   deviance  conditionabacus  0.4142401 0.13779835  3.006132
## 16    2      linr2      value.scale  0.1916233 0.06961279  2.752702
## 17    2      linr2   verbalwm.scale  0.1457604 0.07140294  2.041378
## 18    2      linr2     ravens.scale  0.3187100 0.07676868  4.151563
## 19    2      linr2  conditionabacus  0.4339181 0.13400545  3.238063
## 20    2        ans   verbalwm.scale  0.1740391 0.07312008  2.380182
## 21    2        ans     ravens.scale  0.3182655 0.07901458  4.027934
## 22    2        ans  conditionabacus  0.4344335 0.13760369  3.157135
## 23    2 ordinality   verbalwm.scale  0.1641903 0.07171437  2.289503
## 24    2 ordinality     ravens.scale  0.3052984 0.07787643  3.920293
## 25    2 ordinality  conditionabacus  0.4118030 0.13615714  3.024469
## 26    3   deviance mental.rot.scale  0.2653387 0.07480796  3.546931
## 27    3   deviance     ravens.scale  0.2884607 0.07308947  3.946679
## 28    3   deviance  conditionabacus  0.3285138 0.12625532  2.601979
## 29    3      linr2      value.scale  0.1744391 0.06581358  2.650503
## 30    3      linr2 mental.rot.scale  0.2190144 0.07509538  2.916482
## 31    3      linr2  spatialwm.scale  0.1468878 0.06588906  2.229320
## 32    3      linr2     ravens.scale  0.2881506 0.07133822  4.039218
## 33    3      linr2  conditionabacus  0.2783666 0.12466335  2.232947
## 34    3        ans mental.rot.scale  0.2607191 0.07407022  3.519891
## 35    3        ans     ravens.scale  0.2796418 0.07308300  3.826360
## 36    3        ans  conditionabacus  0.3183353 0.12554495  2.535628
## 37    3 ordinality mental.rot.scale  0.2589200 0.07393976  3.501769
## 38    3 ordinality     ravens.scale  0.2807431 0.07272759  3.860200
## 39    3 ordinality  conditionabacus  0.3196736 0.12493894  2.558639
##         p.value
## 1  2.065866e-02
## 2  2.522071e-02
## 3  6.503401e-03
## 4  4.182095e-02
## 5  1.029879e-02
## 6  4.848038e-02
## 7  1.072408e-02
## 8  1.703141e-02
## 9  4.411998e-02
## 10 1.886275e-02
## 11 4.069567e-02
## 12 1.197156e-02
## 13 1.835480e-02
## 14 2.008624e-04
## 15 3.057078e-03
## 16 6.568543e-03
## 17 4.279644e-02
## 18 5.267619e-05
## 19 1.453375e-03
## 20 1.844511e-02
## 21 8.567010e-05
## 22 1.894650e-03
## 23 2.330961e-02
## 24 1.291290e-04
## 25 2.886892e-03
## 26 5.030810e-04
## 27 1.155482e-04
## 28 1.008181e-02
## 29 8.791686e-03
## 30 4.015162e-03
## 31 2.709501e-02
## 32 8.086261e-05
## 33 2.685006e-02
## 34 5.541413e-04
## 35 1.824602e-04
## 36 1.212612e-02
## 37 5.896754e-04
## 38 1.604113e-04
## 39 1.137530e-02
pca.models$term <- factor(pca.models$term, 
                          levels = c("(Intercept)",
                                     "value.scale", "mental.rot.scale",
                                     "spatialwm.scale", "verbalwm.scale",
                                     "ravens.scale","age", "conditionabacus"), 
                          labels = c("Intercept", 
                                     "Predictor", "Mental Rotation",
                                     "Spatial WM", "Verbal WM",
                                     "Raven's", "Age", "Intervention"))

pca.models$measure <- factor(pca.models$measure,
                              levels = c("ans","deviance","linr2","ordinality"),
                              labels = c("ANS","PAE","Linear r^2", 
                                         "Ordinality"))
  
qplot(term, estimate, 
      fill = term, 
      ymin = estimate - std.error, ymax = estimate + std.error,
      geom = c("bar", "linerange"), stat = "identity", 
      facets = measure ~ year, 
      data=filter(pca.models, term != "Intercept")) + 
  geom_hline(yintercept=0, lty=2) + 
  xlab("Predictor") + 
  ylab("Standardized Beta Weight") + 
  scale_fill_discrete(name="Predictor") +
  theme(axis.text.x = element_text(angle = 90, vjust=.5, hjust = 1)) 

Additional code to double check for no glitches in above code

Now for the guts of one of the models, just an example from year 0, ANS. This uses ANOVA for model comparison.

## run models
# note, the filtering expression gets the model out of the same data frame that we used for the figure
# std.models <- md %>%
#   group_by(year, measure) %>%
#   filter(!(year == 0 & measure != "ans")) %>%
#   do(tidy(lm(standardized.scale ~ value.scale + 
#                mental.rot.scale + 
#                verbalwm.scale + 
#                spatialwm.scale + 
#                ravens.scale + 
#                age + condition, data = .)))

model1 <- lm(standardized.scale ~ value.scale + mental.rot.scale + verbalwm.scale + 
              spatialwm.scale + ravens.scale +  age + condition, 
             data = filter(md, year == 0, measure == "ans", 
                           complete.cases(value.scale)))
model2 <- lm(standardized.scale ~ mental.rot.scale + verbalwm.scale + 
              spatialwm.scale + ravens.scale +  age + condition, 
             data = filter(md, year == 0, measure == "ans", 
                           complete.cases(value.scale)))
summary(model1)
## 
## Call:
## lm(formula = standardized.scale ~ value.scale + mental.rot.scale + 
##     verbalwm.scale + spatialwm.scale + ravens.scale + age + condition, 
##     data = filter(md, year == 0, measure == "ans", complete.cases(value.scale)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.13977 -0.63214 -0.01961  0.62713  2.29587 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      -0.785654   0.979309  -0.802  0.42380   
## value.scale      -0.238260   0.078696  -3.028  0.00295 **
## mental.rot.scale -0.150685   0.082798  -1.820  0.07095 . 
## verbalwm.scale    0.182732   0.092104   1.984  0.04926 * 
## spatialwm.scale   0.229820   0.094789   2.425  0.01663 * 
## ravens.scale     -0.003256   0.082336  -0.040  0.96852   
## age               0.126252   0.146292   0.863  0.38964   
## conditionabacus   0.030503   0.151808   0.201  0.84105   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9024 on 137 degrees of freedom
##   (17 observations deleted due to missingness)
## Multiple R-squared:  0.2187, Adjusted R-squared:  0.1788 
## F-statistic: 5.479 on 7 and 137 DF,  p-value: 1.438e-05
summary(model2)
## 
## Call:
## lm(formula = standardized.scale ~ mental.rot.scale + verbalwm.scale + 
##     spatialwm.scale + ravens.scale + age + condition, data = filter(md, 
##     year == 0, measure == "ans", complete.cases(value.scale)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.40108 -0.61155 -0.01836  0.60736  2.33124 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      -1.178042   0.999003  -1.179  0.24034   
## mental.rot.scale -0.140436   0.085141  -1.649  0.10133   
## verbalwm.scale    0.156911   0.094382   1.663  0.09868 . 
## spatialwm.scale   0.288327   0.095505   3.019  0.00302 **
## ravens.scale      0.009648   0.084623   0.114  0.90939   
## age               0.182500   0.149340   1.222  0.22377   
## conditionabacus   0.069035   0.155685   0.443  0.65815   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9287 on 138 degrees of freedom
##   (17 observations deleted due to missingness)
## Multiple R-squared:  0.1665, Adjusted R-squared:  0.1302 
## F-statistic: 4.593 on 6 and 138 DF,  p-value: 0.0002778
anova(model1, model2)
## Analysis of Variance Table
## 
## Model 1: standardized.scale ~ value.scale + mental.rot.scale + verbalwm.scale + 
##     spatialwm.scale + ravens.scale + age + condition
## Model 2: standardized.scale ~ mental.rot.scale + verbalwm.scale + spatialwm.scale + 
##     ravens.scale + age + condition
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1    137 111.55                                
## 2    138 119.01 -1   -7.4635 9.1663 0.002946 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# also note that p value for the anova is the same as the coefficient p value in the models data frame
filter(std.models, year == 0, measure == "ANS")
## Source: local data frame [8 x 7]
## Groups: year, measure
## 
##   year measure            term     estimate  std.error  statistic
## 1    0     ANS       Intercept -0.785653894 0.97930903 -0.8022533
## 2    0     ANS       Predictor -0.238260490 0.07869646 -3.0275885
## 3    0     ANS Mental Rotation -0.150684986 0.08279761 -1.8199195
## 4    0     ANS       Verbal WM  0.182732396 0.09210359  1.9839879
## 5    0     ANS      Spatial WM  0.229819917 0.09478945  2.4245306
## 6    0     ANS         Raven's -0.003255631 0.08233557 -0.0395410
## 7    0     ANS             Age  0.126251692 0.14629239  0.8630093
## 8    0     ANS    Intervention  0.030503172 0.15180823  0.2009323
##       p.value
## 1 0.423795896
## 2 0.002946144
## 3 0.070954268
## 4 0.049255685
## 5 0.016631616
## 6 0.968516632
## 7 0.389640626
## 8 0.841049533

Year 0 Analyses by JS, to check betas

## 
## Call:
## lm(formula = standardized.scale ~ ans.scale + mental.rot.scale + 
##     verbalwm.scale + spatialwm.scale + ravens.scale + age + condition, 
##     data = yc0)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.15507 -0.63666 -0.01975  0.63162  2.31229 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      -0.858627   0.987588  -0.869  0.38614   
## ans.scale        -0.239697   0.079171  -3.028  0.00295 **
## mental.rot.scale -0.146649   0.080580  -1.820  0.07095 . 
## verbalwm.scale    0.178904   0.090174   1.984  0.04926 * 
## spatialwm.scale   0.231096   0.095316   2.425  0.01663 * 
## ravens.scale     -0.003261   0.082482  -0.040  0.96852   
## age               0.127154   0.147338   0.863  0.38964   
## conditionabacus   0.030721   0.152894   0.201  0.84105   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9088 on 137 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.2187, Adjusted R-squared:  0.1788 
## F-statistic: 5.479 on 7 and 137 DF,  p-value: 1.438e-05
## 
## Call:
## lm(formula = standardized.scale ~ mental.rot.scale + verbalwm.scale + 
##     spatialwm.scale + ravens.scale + age + condition, data = yc0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4182 -0.6159 -0.0185  0.6117  2.3479 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      -1.257841   1.007289  -1.249  0.21387   
## mental.rot.scale -0.136675   0.082861  -1.649  0.10133   
## verbalwm.scale    0.153623   0.092405   1.663  0.09868 . 
## spatialwm.scale   0.289929   0.096035   3.019  0.00302 **
## ravens.scale      0.009665   0.084773   0.114  0.90939   
## age               0.183804   0.150407   1.222  0.22377   
## conditionabacus   0.069529   0.156799   0.443  0.65815   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9353 on 138 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.1665, Adjusted R-squared:  0.1302 
## F-statistic: 4.593 on 6 and 138 DF,  p-value: 0.0002778
## Analysis of Variance Table
## 
## Model 1: standardized.scale ~ ans.scale + mental.rot.scale + verbalwm.scale + 
##     spatialwm.scale + ravens.scale + age + condition
## Model 2: standardized.scale ~ mental.rot.scale + verbalwm.scale + spatialwm.scale + 
##     ravens.scale + age + condition
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1    137 113.15                                
## 2    138 120.72 -1   -7.5706 9.1663 0.002946 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Longitudial growth models

Let’s start by plotting some growth curves for standardized math measures.

qplot(year, standardized, facets=~subnum, 
      geom="line",
      data=subset(d, subnum %in% unique(subnum)[1:30])) + 
  ylim(c(0,1)) + 
  ylab("Standardized Math Tests (WIAT & WJ)") + 
  xlab("Year")

Now add some models. Check whether linear or quadratic is more useful.

preds <- c("standardized","year","subnum")
d.cc <- d[complete.cases(d[,preds]),preds]

mod1 <- lmer(standardized ~ year + (year | subnum), 
            data=d.cc)
mod2 <- lmer(standardized ~ year + I(year^2) + (year + I(year^2) | subnum), 
            data=d.cc)
anova(mod1,mod2)
## refitting model(s) with ML (instead of REML)
## Data: d.cc
## Models:
## mod1: standardized ~ year + (year | subnum)
## mod2: standardized ~ year + I(year^2) + (year + I(year^2) | subnum)
##      Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## mod1  6 -1931.1 -1903.4 971.54  -1943.1                         
## mod2 10 -1928.2 -1882.0 974.10  -1948.2 5.1145      4     0.2758

No gain for this DV (standardized tests) in adding a quadratic term. Let’s drop it then.

Here’s a basic model showing individaul variability plotted by the same model for all participants. It’s a basic linear model and it doesn’t do that badly.

The approach we’re going to follow is to try and see gains by adding predictors to this model to better capture individual variation. (Note that if we use individualized growth terms in a mixed model we already fit pretty much all the variance so there’s not much else we can do, so we are using regular linear models here).

preds <- c("standardized","year","subnum")
d.cc <- d[complete.cases(d[,preds]),preds]

mod <- lm(standardized ~ year, 
            data=d.cc)
d.cc$mod <- predict(mod)

qplot(year, standardized, facets=~subnum, 
      geom="point",
      data=subset(d.cc, subnum %in% unique(subnum)[1:30])) + 
  geom_line(aes(x=year, y=mod),lty=2) + 
  ylim(c(0,1)) + 
  ylab("Standardized Math Tests (WIAT & WJ)") + 
  xlab("Year")

Let’s try adding concurrent ANS as a predictor. As you can see, it does very little here if we have a single coefficient.

(Note that there are some missing values for ANS from year 0, so that’s why curves are truncated).

preds <- c("standardized","year","subnum","ans")
d.cc <- d[complete.cases(d[,preds]),preds]

mod <- lm(standardized ~ year,
            data=d.cc)
mod2 <- lm(standardized ~ year + ans, 
            data=d.cc)
d.cc$mod <- predict(mod)
d.cc$mod2 <- predict(mod2)
anova(mod,mod2)
## Analysis of Variance Table
## 
## Model 1: standardized ~ year
## Model 2: standardized ~ year + ans
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    710 5.2014                                  
## 2    709 5.0825  1   0.11887 16.582 5.184e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qplot(year, standardized, facets=~subnum, 
      geom="point",
      data=subset(d.cc, subnum %in% unique(subnum)[1:30])) + 
  geom_line(aes(x=year, y=mod),lty=2) + 
  geom_line(aes(x=year, y=mod2),lty=1,col="red") + 
  ylim(c(0,1)) + 
  ylab("Standardized Math Tests (WIAT & WJ)") + 
  xlab("Year")

Multiple coefficients, one for each year, and you start to see some teensy-tiny shape differences between curves. (I added abacus condition here as well).

The dashed line is without ANS, the solid is with.

preds <- c("standardized","year","subnum","condition","ans")
d.cc <- d[complete.cases(d[,preds]),preds]

mod <- lm(standardized ~ factor(year) * condition,
            data=d.cc)
mod2 <- lm(standardized ~ factor(year) * condition + 
             factor(year) * ans - ans, 
            data=d.cc)
d.cc$mod <- predict(mod)
d.cc$mod2 <- predict(mod2)
anova(mod,mod2)
## Analysis of Variance Table
## 
## Model 1: standardized ~ factor(year) * condition
## Model 2: standardized ~ factor(year) * condition + factor(year) * ans - 
##     ans
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    704 5.1389                                  
## 2    700 4.9528  4   0.18608 6.5748 3.386e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qplot(year, standardized, facets=~subnum, 
      geom="point",col=condition,
      data=subset(d.cc, subnum %in% unique(subnum)[1:30])) + 
  geom_line(aes(x=year, y=mod),lty=2) + 
  geom_line(aes(x=year, y=mod2),lty=1) + 
  ylim(c(0,1)) + 
  ylab("Standardized Math Tests (WIAT & WJ)") + 
  xlab("Year")

Just for kicks, let’s add everything into the mix and see how we do.

We are doing a bunch better in making customized predictions. In terms of significant effects, we get:

  • arithmetic (big predictor of - you guessed it - standardized arithmetic)
  • abacus condition (whew)
  • mental rotation
  • raven’s
  • verbal wm (but not spatial)

Nothing for ANS. Note that these predictions get noticeably worse when you leave out arithmetic scores. Domain knowledge is a big predictor.

preds <- c("standardized","year","subnum","condition","ans","mental.rot","spatialwm","verbalwm","ravens","age")
d.cc <- d[complete.cases(d[,preds]),preds]

mod <- lm(standardized ~ factor(year) * condition,
            data=d.cc)
mod2 <- lm(standardized ~ factor(year) + condition + 
             age + mental.rot + 
             spatialwm + ravens + verbalwm + ans, 
            data=d.cc)
d.cc$mod <- predict(mod)
d.cc$mod2 <- predict(mod2)

summary(mod2)
## 
## Call:
## lm(formula = standardized ~ factor(year) + condition + age + 
##     mental.rot + spatialwm + ravens + verbalwm + ans, data = d.cc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32374 -0.05260 -0.00326  0.05339  0.20364 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.076577   0.045052   1.700 0.089649 .  
## factor(year)1    0.078369   0.010729   7.305 8.01e-13 ***
## factor(year)2    0.177043   0.011678  15.161  < 2e-16 ***
## factor(year)3    0.280678   0.012651  22.185  < 2e-16 ***
## conditionabacus  0.013745   0.006179   2.224 0.026461 *  
## age              0.004044   0.005990   0.675 0.499858    
## mental.rot       0.069274   0.019269   3.595 0.000348 ***
## spatialwm        0.005252   0.003459   1.518 0.129385    
## ravens           0.084440   0.016781   5.032 6.27e-07 ***
## verbalwm         0.013164   0.003741   3.519 0.000463 ***
## ans             -0.078625   0.031119  -2.527 0.011748 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07883 on 662 degrees of freedom
## Multiple R-squared:  0.7316, Adjusted R-squared:  0.7276 
## F-statistic: 180.5 on 10 and 662 DF,  p-value: < 2.2e-16
qplot(year, standardized, facets=~subnum, 
      geom="point",col=condition,
      data=subset(d.cc, subnum %in% unique(subnum)[1:30])) + 
  geom_line(aes(x=year, y=mod),lty=2) + 
  geom_line(aes(x=year, y=mod2),lty=1) + 
  ylim(c(0,1)) + 
  ylab("Standardized Math Tests (WIAT & WJ)") + 
  xlab("Year")

A final analysis: Let’s try this by lagging the predictors. We use the year before to predict the current year - so this is trying to use e.g. year 0’s ANS to see if it predicts growth over standard predictions in year 1.

Note that I added gender to the baseline predictor set so we will see only lag differences.

Define a function to do the lagging.

lag <- function(x, p) {
  for (y in x$year) {
    if (y == 0) {
      eval(parse(text=paste("x$",p,".lag[x$year==0] <- 0",sep="")))
      } 
    
    if (y != 0 & ((y-1) %in% x$year)) {
      eval(parse(text=paste("x$",p,".lag[x$year==y] <- x$",
                            p,"[x$year==y-1]",sep="")))
      } else {
      eval(parse(text=paste("x$",p,".lag[x$year==y] <- NA",sep="")))      
      }
    }
  return(x)
  }

Now do the actual predictions. We are actually missing too many observations to really do this right without some serious imputation. We lose around a third of our data, as you can see from the plot.

Because of the missing data problem, do this with just ANS, to see if anything else is different. Results:

  • ANS is a predictor
  • lagged ANS isn’t
  • None of this is visible in the plot because the differences are too small, so I left out the plot.